"""
Module containing the functions to handle the Moe&Distrefano data
This class object is an extension to the population grid object
TODO: The values from M&S are from bincenters, but our method makes it that the bincenter value is used as an edge value. So e.g. taking an eccentricity of 0.04 wll give 0 because its lower than the 0.05 (minimum) point. But it should probably give a value, even if we choose eccentricity = 0
"""
# pylint: disable=E1101
import copy
import gc
import json
import os
import py_rinterpolate
from binarycpython.utils import moe_di_stefano_2017_data
from binarycpython.utils.dicts import normalize_dict, update_dicts
from binarycpython.utils.population_extensions.distribution_functions import (
LOG_LN_CONVERTER,
Moecache,
)
[docs]class Moe_di_Stefano_2017:
"""
Extension to the population grid object that contains functionality to handle handle the Moe & distefano distributions
"""
def __init__(self, **kwargs):
"""
Init function for the Moe_di_Stefano_2017 class
"""
self.moe_distefano_options_defaults_dict = {
"JSON": {
"value": None,
"description": "Parameter that stores the interpolation data",
},
"resolutions": {
"value": {
"M": [
20, # M1
20, # M2 (i.e. q)
0, # M3 currently unused
0, # M4 currently unused
],
"logP": [
20, # P2 (binary period)
0, # P3 (triple period) currently unused
0, # P4 (quadruple period) currently unused
],
"ecc": [
10, # e (binary eccentricity)
0, # e2 (triple eccentricity) currently unused
0, # e3 (quadruple eccentricity) currently unused
],
},
"description": "Dictionary that stores the resolutions for each of the variables that are used to sample the systems. The structure is as follows: resolutions = {'M': [M1 res, M2 res, ..], 'logP': [logP1 res, logP2 res, ...], 'ecc': [ecc1 res, ecc2 res, ...]}",
},
"samplerfuncs": {
"value": {
"M": [None, None, None, None],
"logP": [None, None, None],
"ecc": [None, None, None],
},
"description": "sampler functions to each of the parameters. NEEDS UPDATING",
},
"IMF_distribution": {
"value": "kroupa2001",
"description": "IMF choice for the M&S sampling. Currently only supporting 'kroupa2001': Kroupa 2001 IMF and 'chabrier2003': Chabrier 2003 IMF.",
},
"ranges": {
"value": {
# stellar masses (Msun)
"M": [
self.minimum_stellar_mass()
* 1.05, # 0.08 is a tad bit above the minimum mass. Don't sample at 0.07, otherwise the first row of q values will have a phasevol of 0. Anything higher is fine.
80.0, # (rather arbitrary) upper mass cutoff
],
"q": [
None, # artificial qmin : set to None to use default
None, # artificial qmax : set to None to use default
],
"logP": [0.0, 8.0], # 0 = log10(1 day) # 8 = log10(10^8 days)
"ecc": [0.0, 0.99],
},
"description": "Ranges for the parameters that are sampled through the M&S distributions. Input is as follows: resolutions = {'M': [M1 lower bound, M2 upper bound], 'q': [q lower bound, q upper bound]}",
},
"Mmin": {
"value": self.minimum_stellar_mass(),
"description": "Minimum stellar mass that is regarded a star",
},
"multiplicity_model": {
"value": "Poisson",
"description": """
multiplicity model (as a function of log10M1)
You can use 'Poisson' which uses the system multiplicity
given by Moe and maps this to single/binary/triple/quad
fractions.
Alternatively, 'data' takes the fractions directly
from the data, but then triples and quadruples are
combined (and there are NO quadruples).
""",
},
"multiplicity_modulator": {
"value": [
1, # single
1, # binary
0, # triple
0, # quadruple
],
"description": """
[single, binary, triple, quadruple]
e.g. [1,0,0,0] for single stars only
[0,1,0,0] for binary stars only
defaults to [1,1,0,0] i.e. singles and binaries
""",
},
"normalize_multiplicities": {
"value": "merge",
"description": """
'norm': normalise so the whole population is 1.0
after implementing the appropriate fractions
S/(S+B+T+Q), B/(S+B+T+Q), T/(S+B+T+Q), Q/(S+B+T+Q)
given a mix of multiplicities, you can either (noting that
here (S,B,T,Q) = appropriate modulator * model(S,B,T,Q) )
note: if you only set one multiplicity_modulator
to 1, and all the others to 0, then normalising
will mean that you effectively have the same number
of stars as single, binary, triple or quad (whichever
is non-zero) i.e. the multiplicity fraction is ignored.
This is probably not useful except for
testing purposes or comparing to old grids.
'raw' : stick to what is predicted, i.e.
S/(S+B+T+Q), B/(S+B+T+Q), T/(S+B+T+Q), Q/(S+B+T+Q)
without normalisation
(in which case the total probability < 1.0 unless
all you use single, binary, triple and quadruple)
'merge' : e.g. if you only have single and binary,
add the triples and quadruples to the binaries, so
binaries represent all multiple systems
...
*** this is canonical binary population synthesis ***
It only takes the maximum multiplicity into account,
i.e. it doesn't multiply the resulting array by the multiplicity modulator again.
This prevents the resulting array to always be 1 if only 1 multiplicity modulator element is nonzero
Note: if multiplicity_modulator == [1,1,1,1]. this option does nothing (equivalent to 'raw').
""",
},
"q_low_extrapolation_method": {
"value": "flat",
"description": """
q extrapolation (below 0.15) method
none
flat
linear2
plaw2
nolowq
""",
},
"q_high_extrapolation_method": {
"value": "flat",
"description": "Same as q_low_extrapolation_method",
},
}
[docs] def set_moe_di_stefano_settings(self, options=None):
"""
Function to set user input configurations for the Moe & di Stefano methods
If nothing is passed then we just use the default options
"""
if not options:
options = {}
# Take the option dictionary that was given and override.
options = update_dicts(self.population_options["Moe2017_options"], options)
self.population_options["Moe2017_options"] = copy.deepcopy(options)
# Write options to a file
os.makedirs(
os.path.join(self.population_options["tmp_dir"], "moe_distefano"),
exist_ok=True,
)
with open(
os.path.join(
os.path.join(self.population_options["tmp_dir"], "moe_distefano"),
"moeopts.dat",
),
"w",
encoding="utf-8",
) as f:
f.write(
json.dumps(
self.population_options["Moe2017_options"],
indent=4,
ensure_ascii=False,
)
)
f.close()
def _load_moe_di_stefano_data(self):
"""
Function to load the moe & di stefano data
"""
# Only if the grid is loaded and Moecache contains information
if not self.population_options["_loaded_Moe2017_data"]: # and not Moecache:
if self.population_options["_Moe2017_JSON_data"]:
# Use the existing (perhaps modified) JSON data
json_data = self.population_options["_Moe2017_JSON_data"]
else:
# Load the JSON data from a file
json_data = self.get_moe_di_stefano_dataset(
self.population_options["Moe2017_options"],
verbosity=self.population_options["verbosity"],
)
# entry of log10M1 is a list containing 1 dict.
# We can take the dict out of the list
if isinstance(json_data["log10M1"], list):
json_data["log10M1"] = json_data["log10M1"][0]
# save this data in case we want to modify it later
self.population_options["_Moe2017_JSON_data"] = json_data
# Get all the masses
logmasses = sorted(json_data["log10M1"].keys())
if not logmasses:
msg = "The table does not contain masses."
self.vb_debug(
"\tMoe_di_Stefano_2017: {}".format(msg),
)
raise ValueError(msg)
# Write to file
os.makedirs(
os.path.join(self.population_options["tmp_dir"], "moe_distefano"),
exist_ok=True,
)
with open(
os.path.join(
os.path.join(self.population_options["tmp_dir"], "moe_distefano"),
"moe.log",
),
"w",
encoding="utf-8",
) as logfile:
logfile.write("log₁₀Masses(M☉) {}\n".format(logmasses))
# Get all the periods and see if they are all consistently present
logperiods = []
for logmass in logmasses:
if not logperiods:
logperiods = sorted(json_data["log10M1"][logmass]["logP"].keys())
dlog10P = float(logperiods[1]) - float(logperiods[0])
current_logperiods = sorted(json_data["log10M1"][logmass]["logP"])
if logperiods != current_logperiods:
msg = (
"Period values are not consistent throughout the dataset logperiods = "
+ " ".join(str(x) for x in logperiods)
+ "\nCurrent periods = "
+ " ".join(str(x) for x in current_logperiods)
)
self.vb_debug(
"\tMoe_di_Stefano_2017: {}".format(msg),
)
raise ValueError(msg)
############################################################
# log10period binwidth : of course this assumes a fixed
# binwidth, so we check for this too.
for i in range(len(current_logperiods) - 1):
if not dlog10P == (
float(current_logperiods[i + 1]) - float(current_logperiods[i])
):
msg = "Period spacing is not consistent throughout the dataset"
self.vb_debug(
"\tMoe_di_Stefano_2017: {}".format(msg),
)
raise ValueError(msg)
# save the logperiods list in the cache:
# this is used in the renormalization integration
Moecache["logperiods"] = logperiods
# Write to file
os.makedirs(
os.path.join(self.population_options["tmp_dir"], "moe_distefano"),
exist_ok=True,
)
with open(
os.path.join(
self.population_options["tmp_dir"], "moe_distefano", "moe.log"
),
"a",
encoding="utf-8",
) as logfile:
logfile.write("log₁₀Periods(days) {}\n".format(logperiods))
# Fill the global dict
for logmass in logmasses:
# Create the multiplicity table
if not Moecache.get("multiplicity_table", None):
Moecache["multiplicity_table"] = []
# multiplicity as a function of primary mass
Moecache["multiplicity_table"].append(
[
float(logmass),
json_data["log10M1"][logmass]["f_multi"],
json_data["log10M1"][logmass]["single star fraction"],
json_data["log10M1"][logmass]["binary star fraction"],
json_data["log10M1"][logmass]["triple/quad star fraction"],
]
)
############################################################
# a small log10period which we can shift just outside the
# table to force integration out there to zero
epslog10P = 1e-8 * dlog10P
############################################################
# loop over either binary or triple-outer periods
first = 1
# Go over the periods
for logperiod in logperiods:
############################################################
# distributions of binary and triple star fractions
# as a function of mass, period.
#
# Note: these should be per unit log10P, hence we
# divide by dlog10P
if first:
first = 0
# Create the multiplicity table
if not Moecache.get("period_distributions", None):
Moecache["period_distributions"] = []
############################################################
# lower bound the period distributions to zero probability
Moecache["period_distributions"].append(
[
float(logmass),
float(logperiod) - 0.5 * dlog10P - epslog10P,
0.0,
0.0,
]
)
Moecache["period_distributions"].append(
[
float(logmass),
float(logperiod) - 0.5 * dlog10P,
json_data["log10M1"][logmass]["logP"][logperiod][
"normed_bin_frac_p_dist"
]
/ dlog10P,
json_data["log10M1"][logmass]["logP"][logperiod][
"normed_tripquad_frac_p_dist"
]
/ dlog10P,
]
)
Moecache["period_distributions"].append(
[
float(logmass),
float(logperiod),
json_data["log10M1"][logmass]["logP"][logperiod][
"normed_bin_frac_p_dist"
]
/ dlog10P,
json_data["log10M1"][logmass]["logP"][logperiod][
"normed_tripquad_frac_p_dist"
]
/ dlog10P,
]
)
############################################################
# distributions as a function of mass, period, q
#
# First, get a list of the qs given by Moe
#
qs = sorted(json_data["log10M1"][logmass]["logP"][logperiod]["q"])
# Fill the data and 'normalise'
qdata = self.fill_data(
qs, json_data["log10M1"][logmass]["logP"][logperiod]["q"]
)
# Create the multiplicity table
if not Moecache.get("q_distributions", None):
Moecache["q_distributions"] = []
for q in qs:
Moecache["q_distributions"].append(
[float(logmass), float(logperiod), float(q), qdata[q]]
)
############################################################
# eccentricity distributions as a function of mass, period, ecc
eccs = sorted(json_data["log10M1"][logmass]["logP"][logperiod]["e"])
# Fill the data and 'normalise'
ecc_data = self.fill_data(
eccs, json_data["log10M1"][logmass]["logP"][logperiod]["e"]
)
# Create the multiplicity table
if not Moecache.get("ecc_distributions", None):
Moecache["ecc_distributions"] = []
for ecc in eccs:
Moecache["ecc_distributions"].append(
[
float(logmass),
float(logperiod),
float(ecc),
ecc_data[ecc],
]
)
############################################################
# upper bound the period distributions to zero probability
Moecache["period_distributions"].append(
[
float(logmass),
float(logperiods[-1])
+ 0.5 * dlog10P, # TODO: why this shift? to center it?
json_data["log10M1"][logmass]["logP"][logperiods[-1]][
"normed_bin_frac_p_dist"
]
/ dlog10P,
json_data["log10M1"][logmass]["logP"][logperiods[-1]][
"normed_tripquad_frac_p_dist"
]
/ dlog10P,
]
)
Moecache["period_distributions"].append(
[
float(logmass),
float(logperiods[-1]) + 0.5 * dlog10P + epslog10P,
0.0,
0.0,
]
)
self.vb_debug(
"\tMoe_di_Stefano_2017: Length period_distributions table: {}".format(
len(Moecache["period_distributions"])
),
)
self.vb_debug(
"\tMoe_di_Stefano_2017: Length multiplicity table: {}".format(
len(Moecache["multiplicity_table"])
),
)
self.vb_debug(
"\tMoe_di_Stefano_2017: Length q table: {}".format(
len(Moecache["q_distributions"])
),
)
self.vb_debug(
"\tMoe_di_Stefano_2017: Length ecc table: {}".format(
len(Moecache["ecc_distributions"])
),
)
# Write to log file
os.makedirs(
os.path.join(self.population_options["tmp_dir"], "moe_distefano"),
exist_ok=True,
)
with open(
os.path.join(
os.path.join(self.population_options["tmp_dir"], "moe_distefano"),
"moecache.json",
),
"w",
encoding="utf-8",
) as cache_filehandle:
cache_filehandle.write(
json.dumps(Moecache, indent=4, ensure_ascii=False)
)
# Signal that the data has been loaded
self.population_options["_loaded_Moe2017_data"] = True
def _set_moe_di_stefano_distributions(self):
"""
Function to set the Moe & di Stefano distribution
"""
##############
# Handle multiplicity loop if we are doing grid sampling
if not self.population_options["evolution_type"] == "monte_carlo":
############################################################
# first, the multiplicity, this is 1,2,3,4, ...
# for singles, binaries, triples, quadruples, ...
max_multiplicity = self.get_max_multiplicity(
self.population_options["Moe2017_options"]["multiplicity_modulator"]
)
self.vb_debug(
"\tMoe_di_Stefano_2017: Max multiplicity = {}".format(max_multiplicity),
)
# Multiplicity
self.add_sampling_variable(
name="multiplicity",
parameter_name="multiplicity",
longname="multiplicity",
valuerange=[1, max_multiplicity],
samplerfunc="self.const_int(1, {n}, {n})".format(n=max_multiplicity),
precode='self.population_options["multiplicity"] = multiplicity; self.bse_options["multiplicity"] = multiplicity; options={}'.format(
self.population_options["Moe2017_options"]
),
condition="({}[int(multiplicity)-1] > 0)".format(
str(
self.population_options["Moe2017_options"][
"multiplicity_modulator"
]
)
),
gridtype="discrete",
probdist=1,
)
####################################
# Set up the sampling variables
################
# Primary mass: M_1
self.add_sampling_variable(
name="lnM_1",
parameter_name="M_1",
longname="Primary mass",
samplerfunc=self.population_options["Moe2017_options"]["samplerfuncs"]["M"][
0
]
or "self.const_linear(np.log({}), np.log({}), {})".format(
self.population_options["Moe2017_options"]["ranges"]["M"][0],
self.population_options["Moe2017_options"]["ranges"]["M"][1],
self.population_options["Moe2017_options"]["resolutions"]["M"][0],
),
valuerange=[
"np.log({})".format(
self.population_options["Moe2017_options"]["ranges"]["M"][0]
),
"np.log({})".format(
self.population_options["Moe2017_options"]["ranges"]["M"][1]
),
],
gridtype="centred",
dphasevol="dlnM_1",
precode='M_1 = np.exp(lnM_1); options["M_1"]=M_1',
probdist="self.Moe_di_Stefano_2017_pdf({{{}, {}, {}}}, verbosity=self.population_options['verbosity'])['total_probdens'] if multiplicity == 1 else 1".format(
str(dict(self.population_options["Moe2017_options"]))[1:-1],
"'multiplicity': multiplicity",
"'M_1': M_1",
),
###########
# Monte-Carlo sampling related
# To indicate that this is a branchpoint and we want to execute this when the multiplicity==1
branchpoint=1,
branchcode="multiplicity == 1",
)
################
# Handle binaries
if max_multiplicity >= 2:
################
# Orbital period inner binary: orbital_period
self.add_sampling_variable(
name="log10per",
parameter_name="orbital_period",
longname="log10(Orbital_Period)",
probdist=1.0,
condition='(self.population_options["multiplicity"] >= 2)',
gridtype="centred",
dphasevol="({} * dlog10per)".format(LOG_LN_CONVERTER),
valuerange=[
self.population_options["Moe2017_options"]["ranges"]["logP"][0],
self.population_options["Moe2017_options"]["ranges"]["logP"][1],
],
samplerfunc=self.population_options["Moe2017_options"]["samplerfuncs"][
"logP"
][0]
or "self.const_linear({}, {}, {})".format(
self.population_options["Moe2017_options"]["ranges"]["logP"][0],
self.population_options["Moe2017_options"]["ranges"]["logP"][1],
self.population_options["Moe2017_options"]["resolutions"]["logP"][
0
],
),
precode="""orbital_period = 10.0**log10per
qmin={}/M_1
qmax=maximum_mass_ratio_for_RLOF(M_1, orbital_period)
""".format(
self.population_options["Moe2017_options"]["Mmin"]
),
) # TODO: change the maximum_mass_ratio_for_RLOF
################
# Mass ratio inner binary: q / M_2
self.add_sampling_variable(
name="q",
parameter_name="M_2",
longname="Mass ratio",
valuerange=[
self.population_options["Moe2017_options"]["ranges"]["q"][0]
if self.population_options["Moe2017_options"]
.get("ranges", {})
.get("q", None)
else "options['Mmin']/M_1",
self.population_options["Moe2017_options"]["ranges"]["q"][1]
if self.population_options["Moe2017_options"]
.get("ranges", {})
.get("q", None)
else "qmax",
],
probdist=1,
gridtype="centred",
dphasevol="dq",
precode="""
M_2 = q * M_1
sep = calc_sep_from_period(M_1, M_2, orbital_period)
""",
samplerfunc=self.population_options["Moe2017_options"]["samplerfuncs"][
"M"
][1]
or "self.const_linear({}, {}, {})".format(
self.population_options["Moe2017_options"]["ranges"]["q"][0]
if self.population_options["Moe2017_options"]
.get("ranges", {})
.get("q", [None, None])[0]
else "{}/M_1".format(
self.population_options["Moe2017_options"]["Mmin"]
),
self.population_options["Moe2017_options"]["ranges"]["q"][1]
if self.population_options["Moe2017_options"]
.get("ranges", {})
.get("q", [None, None])[1]
else "qmax",
self.population_options["Moe2017_options"]["resolutions"]["M"][1],
),
###########
# Monte-Carlo sampling related
# To indicate that this is a branchpoint and we want to execute this when the multiplicity==2
branchpoint=1
if (
self.population_options["evolution_type"] == "evolution_type"
and self.population_options["Moe2017_options"]["resolutions"][
"ecc"
][0]
== 0
)
else 0,
branchcode="multiplicity == 2"
if (
self.population_options["evolution_type"] == "evolution_type"
and self.population_options["Moe2017_options"]["resolutions"][
"ecc"
][0]
== 0
)
else None,
)
################
# (optional) eccentricity inner binary: ecc
if self.population_options["Moe2017_options"]["resolutions"]["ecc"][0] > 0:
self.add_sampling_variable(
name="ecc",
parameter_name="eccentricity",
longname="Eccentricity",
probdist=1,
gridtype="centred",
dphasevol="decc",
precode="eccentricity=ecc",
valuerange=[
self.population_options["Moe2017_options"]["ranges"]["ecc"][
0
], # Just fail if not defined.
self.population_options["Moe2017_options"]["ranges"]["ecc"][1],
],
samplerfunc=self.population_options["Moe2017_options"][
"samplerfuncs"
]["ecc"][0]
or "self.const_linear({}, {}, {})".format(
self.population_options["Moe2017_options"]["ranges"]["ecc"][
0
], # Just fail if not defined.
self.population_options["Moe2017_options"]["ranges"]["ecc"][1],
self.population_options["Moe2017_options"]["resolutions"][
"ecc"
][0],
),
###########
# Monte-Carlo sampling related
# To indicate that this is a branchpoint and we want to execute this when the multiplicity==2
branchpoint=1
if (self.population_options["evolution_type"] == "evolution_type")
else 0,
branchcode="multiplicity == 2"
if (self.population_options["evolution_type"] == "evolution_type")
else None,
)
################
# Handle triple systems
if max_multiplicity >= 3:
################
# Orbital period hierarchical triple system: orbital_period_triple
self.add_sampling_variable(
name="log10per2",
parameter_name="orbital_period_triple",
longname="log10(Orbital_Period2)",
probdist=1.0,
condition='(self.population_options["multiplicity"] >= 3)',
gridtype="centred",
dphasevol="({} * dlog10per2)".format(LOG_LN_CONVERTER),
valuerange=[
self.population_options["Moe2017_options"]["ranges"]["logP"][0],
self.population_options["Moe2017_options"]["ranges"]["logP"][1],
],
samplerfunc=self.population_options["Moe2017_options"]["samplerfuncs"][
"logP"
][1]
or "self.const_linear({}, {}, {})".format(
self.population_options["Moe2017_options"]["ranges"]["logP"][0],
self.population_options["Moe2017_options"]["ranges"]["logP"][1],
self.population_options["Moe2017_options"]["resolutions"]["logP"][
1
],
),
precode="""orbital_period_triple = 10.0**log10per2
q2min={}/(M_1+M_2)
q2max=maximum_mass_ratio_for_RLOF(M_1+M_2, orbital_period_triple)
""".format(
self.population_options["Moe2017_options"]["Mmin"]
),
)
################
# mass ratio hierarchical triple system: q2 & M_3
# NOTE: the mass ratio is M_outer/M_inner
self.add_sampling_variable(
name="q2",
parameter_name="M_3",
longname="Mass ratio outer/inner",
valuerange=[
self.population_options["Moe2017_options"]["ranges"]["q"][0]
if self.population_options["Moe2017_options"]
.get("ranges", {})
.get("q", None)
else "options['Mmin']/(M_1+M_2)",
self.population_options["Moe2017_options"]["ranges"]["q"][1]
if self.population_options["Moe2017_options"]
.get("ranges", {})
.get("q", None)
else "q2max",
],
probdist=1,
gridtype="centred",
dphasevol="dq2",
precode="""
M_3 = q2 * (M_1 + M_2)
sep2 = calc_sep_from_period((M_1+M_2), M_3, orbital_period_triple)
eccentricity2=0
""",
samplerfunc=self.population_options["Moe2017_options"]["samplerfuncs"][
"M"
][2]
or "self.const_linear({}, {}, {})".format(
self.population_options["Moe2017_options"]["ranges"]["q"][0]
if self.population_options["Moe2017_options"]
.get("ranges", {})
.get("q", None)
else "options['Mmin']/(M_1+M_2)",
self.population_options["Moe2017_options"]["ranges"]["q"][1]
if self.population_options["Moe2017_options"]
.get("ranges", {})
.get("q", None)
else "q2max",
self.population_options["Moe2017_options"]["resolutions"]["M"][2],
),
###########
# Monte-Carlo sampling related
# To indicate that this is a branchpoint and we want to execute this when the multiplicity==3
branchpoint=1
if (
self.population_options["evolution_type"] == "evolution_type"
and self.population_options["Moe2017_options"]["resolutions"][
"ecc"
][1]
== 0
)
else 0,
branchcode="multiplicity == 3"
if (
self.population_options["evolution_type"] == "evolution_type"
and self.population_options["Moe2017_options"]["resolutions"][
"ecc"
][1]
== 0
)
else None,
)
################
# (optional) eccentricity hierarchical triple system: eccentricity_triple
if self.population_options["Moe2017_options"]["resolutions"]["ecc"][1] > 0:
self.add_sampling_variable(
name="ecc2",
parameter_name="eccentricity2",
longname="Eccentricity of the triple",
probdist=1,
gridtype="centred",
dphasevol="decc2",
precode="eccentricity_triple=ecc2",
valuerange=[
self.population_options["Moe2017_options"]["ranges"]["ecc"][
0
], # Just fail if not defined.
self.population_options["Moe2017_options"]["ranges"]["ecc"][1],
],
samplerfunc=self.population_options["Moe2017_options"][
"samplerfuncs"
]["ecc"][1]
or "self.const_linear({}, {}, {})".format(
self.population_options["Moe2017_options"]["ranges"]["ecc"][
0
], # Just fail if not defined.
self.population_options["Moe2017_options"]["ranges"]["ecc"][1],
self.population_options["Moe2017_options"]["resolutions"][
"ecc"
][1],
),
###########
# Monte-Carlo sampling related
# To indicate that this is a branchpoint and we want to execute this when the multiplicity==3
branchpoint=1
if (self.population_options["evolution_type"] == "evolution_type")
else 0,
branchcode="multiplicity == 3"
if (self.population_options["evolution_type"] == "evolution_type")
else None,
)
################
# Handle quadruple systems
if max_multiplicity == 4:
################
# Orbital period hierarchical quadruple system: orbital_period_quadruple
self.add_sampling_variable(
name="log10per3",
parameter_name="orbital_period_quadruple",
longname="log10(Orbital_Period3)",
probdist=1.0,
condition='(self.population_options["multiplicity"] >= 4)',
branchpoint=3
if max_multiplicity > 3
else 0, # Signal here to put a branchpoint if we have a max multiplicity higher than 1.
gridtype="centred",
dphasevol="({} * dlog10per3)".format(LOG_LN_CONVERTER),
valuerange=[
self.population_options["Moe2017_options"]["ranges"]["logP"][0],
self.population_options["Moe2017_options"]["ranges"]["logP"][1],
],
samplerfunc=self.population_options["Moe2017_options"]["samplerfuncs"][
"logP"
][2]
or "self.const_linear({}, {}, {})".format(
self.population_options["Moe2017_options"]["ranges"]["logP"][0],
self.population_options["Moe2017_options"]["ranges"]["logP"][1],
self.population_options["Moe2017_options"]["resolutions"]["logP"][
2
],
),
precode="""orbital_period_quadruple = 10.0**log10per3
q3min={}/(M_3)
q3max=maximum_mass_ratio_for_RLOF(M_3, orbital_period_quadruple)
""".format(
self.population_options["Moe2017_options"]["Mmin"]
),
)
################
# mass ratio hierarchical triple system: q3 & M_4
# NOTE: the mass ratio is M_outer/M_inner
self.add_sampling_variable(
name="q3",
parameter_name="M_4",
longname="Mass ratio outer low/outer high",
valuerange=[
self.population_options["Moe2017_options"]["ranges"]["q"][0]
if self.population_options["Moe2017_options"]
.get("ranges", {})
.get("q", None)
else "options['Mmin']/(M_3)",
self.population_options["Moe2017_options"]["ranges"]["q"][1]
if self.population_options["Moe2017_options"]
.get("ranges", {})
.get("q", None)
else "q3max",
],
probdist=1,
gridtype="centred",
dphasevol="dq3",
precode="""
M_4 = q3 * M_3
sep3 = calc_sep_from_period((M_3), M_4, orbital_period_quadruple)
eccentricity3=0
""",
samplerfunc=self.population_options["Moe2017_options"]["samplerfuncs"][
"M"
][3]
or "self.const_linear({}, {}, {})".format(
self.population_options["Moe2017_options"]["ranges"]["q"][0]
if self.population_options["Moe2017_options"]
.get("ranges", {})
.get("q", None)
else "options['Mmin']/(M_3)",
self.population_options["Moe2017_options"]["ranges"]["q"][1]
if self.population_options["Moe2017_options"]
.get("ranges", {})
.get("q", None)
else "q3max",
self.population_options["Moe2017_options"]["resolutions"]["M"][2],
),
###########
# Monte-Carlo sampling related
# To indicate that this is a branchpoint and we want to execute this when the multiplicity==4
branchpoint=1
if (
self.population_options["evolution_type"] == "evolution_type"
and self.population_options["Moe2017_options"]["resolutions"][
"ecc"
][2]
== 0
)
else 0,
branchcode="multiplicity == 4"
if (
self.population_options["evolution_type"] == "evolution_type"
and self.population_options["Moe2017_options"]["resolutions"][
"ecc"
][2]
== 0
)
else None,
)
################
# (optional) eccentricity hierarchical quadruple system: eccentricity_quadruple
if self.population_options["Moe2017_options"]["resolutions"]["ecc"][2] > 0:
self.add_sampling_variable(
name="ecc3",
parameter_name="eccentricity3",
longname="Eccentricity of the triple+quadruple/outer binary",
probdist=1,
gridtype="centred",
dphasevol="decc3",
precode="eccentricity3=ecc3",
valuerange=[
self.population_options["Moe2017_options"]["ranges"]["ecc"][
0
], # Just fail if not defined.
self.population_options["Moe2017_options"]["ranges"]["ecc"][1],
],
samplerfunc=self.population_options["Moe2017_options"][
"samplerfuncs"
]["ecc"][2]
or "self.const_linear({}, {}, {})".format(
self.population_options["Moe2017_options"]["ranges"]["ecc"][
0
], # Just fail if not defined.
self.population_options["Moe2017_options"]["ranges"]["ecc"][1],
self.population_options["Moe2017_options"]["resolutions"][
"ecc"
][2],
),
###########
# Monte-Carlo sampling related
# To indicate that this is a branchpoint and we want to execute this when the multiplicity==4
branchpoint=1
if (self.population_options["evolution_type"] == "evolution_type")
else 0,
branchcode="multiplicity == 4"
if (self.population_options["evolution_type"] == "evolution_type")
else None,
)
###################
# Wrap up
# Now we are at the last part.
# Here we should combine all the information that we calculate and update the options
# dictionary. This will then be passed to the Moe_di_Stefano_2017_pdf to calculate
# the real probability. The trick we use is to strip the options_dict as a string
# and add some keys to it:
# TODO: we need to add this to each 'final' sampling variable at each metallicity. The stochastic nature of monte-carlo sampling makes that we have changing multipliticies
updated_options = "{{{}, {}, {}, {}, {}, {}, {}, {}, {}, {}, {}, {}}}".format(
str(dict(self.population_options["Moe2017_options"]))[1:-1],
'"multiplicity": multiplicity',
'"M_1": M_1',
'"M_2": M_2',
'"M_3": M_3',
'"M_4": M_4',
'"P": orbital_period',
'"P2": orbital_period_triple',
'"P3": orbital_period_quadruple',
'"ecc": eccentricity',
'"ecc2": eccentricity_triple',
'"ecc3": eccentricity_quadruple',
)
probdist_addition = "self.Moe_di_Stefano_2017_pdf({}, verbosity=self.population_options['verbosity'])['total_probdens']".format(
updated_options
)
# and finally the probability calculator
self.population_options["_sampling_variables"][self._last_sampling_variable()][
"probdist"
] = probdist_addition
self.vb_debug(
"\tMoe_di_Stefano_2017: Added final call to the pdf function",
)
# Signal that the MOE2017 grid has been set
self.population_options["_set_Moe2017_grid"] = True
################################################################################################
[docs] def Moe_di_Stefano_2017(self, options=None):
"""
Function to handle setting the user input settings,
set up the data and load that into interpolators and
then set the distribution functions
Takes a dictionary as its only argument
"""
default_options = {
"apply settings": True,
"setup grid": True,
"load data": True,
"clean cache": False,
"clean load flag": False,
"clean all": False,
}
if not options:
options = {}
options = update_dicts(default_options, options)
# clean cache?
if options["clean all"] or options["clean cache"]:
Moecache.clear()
if options["clean all"] or options["clean load flag"]:
self.population_options["_loaded_Moe2017_data"] = False
# Set the user input
if options["apply settings"]:
self.set_moe_di_stefano_settings(options=options)
# Load the data
if options["load data"]:
self._load_moe_di_stefano_data()
# construct the grid here
if options["setup grid"]:
self._set_moe_di_stefano_distributions()
def _clean_interpolators(self):
"""
Function to clean up the interpolators after a run
We look in the Moecache global variable for items that are interpolators.
Should be called by the general cleanup function AND the thread cleanup function
"""
interpolator_keys = []
for key, value in Moecache.items():
if isinstance(value, py_rinterpolate.Rinterpolate):
interpolator_keys.append(key)
for key in interpolator_keys:
Moecache[key].destroy()
del Moecache[key]
gc.collect()
def _calculate_multiplicity_fraction(self, system_dict):
"""
Function to calculate multiplicity fraction
Makes use of the self.bse_options['multiplicity'] value. If its not set, it will raise an error
population_options['multiplicity_fraction_function'] will be checked for the choice
TODO: add option to put a manual binary fraction in here (solve via negative numbers being the functions)
"""
# Just return 1 if no option has been chosen
if self.population_options["multiplicity_fraction_function"] in [0, "None"]:
self.vb_debug(
"_calculate_multiplicity_fraction: Chosen not to use any multiplicity fraction.",
)
return 1
# Raise an error if the multiplicity is not set
if not system_dict.get("multiplicity", None):
msg = "Multiplicity value has not been set. When using a specific multiplicity fraction function please set the multiplicity"
raise ValueError(msg)
# Go over the chosen options
if self.population_options["multiplicity_fraction_function"] in [
1,
"Arenou2010",
]:
# Arenou 2010 will be used
self.vb_debug(
"_calculate_multiplicity_fraction: Using Arenou 2010 to calculate multiplicity fractions",
)
binary_fraction = self.Arenou2010_binary_fraction(system_dict["M_1"])
multiplicity_fraction_dict = {
1: 1 - binary_fraction,
2: binary_fraction,
3: 0,
4: 0,
}
elif self.population_options["multiplicity_fraction_function"] in [
2,
"Raghavan2010",
]:
# Raghavan 2010 will be used
self.vb_debug(
"_calculate_multiplicity_fraction: Using Raghavan (2010) to calculate multiplicity fractions",
)
binary_fraction = self.raghavan2010_binary_fraction(system_dict["M_1"])
multiplicity_fraction_dict = {
1: 1 - binary_fraction,
2: binary_fraction,
3: 0,
4: 0,
}
elif self.population_options["multiplicity_fraction_function"] in [
3,
"Moe2017",
]:
# We need to check several things now here:
# First, are the options for the MOE2017 grid set? On start it is filled with the default settings
if not self.population_options["Moe2017_options"]:
msg = "The MOE2017 options do not seem to be set properly. The value is {}".format(
self.population_options["Moe2017_options"]
)
raise ValueError(msg)
# Second: is the Moecache filled.
if not Moecache:
self.vb_debug(
"_calculate_multiplicity_fraction: Moecache is empty. It needs to be filled with the data for the interpolators. Loading the data now",
)
# Load the data
self._load_moe_di_stefano_data()
# record the prev value
prev_M1_value_ms = self.population_options["Moe2017_options"].get(
"M_1", None
)
# Set value of M1 of the current system
self.population_options["Moe2017_options"]["M_1"] = system_dict["M_1"]
# Calculate the multiplicity fraction
multiplicity_fraction_list = (
self.Moe_di_Stefano_2017_multiplicity_fractions(
self.population_options["Moe2017_options"],
self.population_options["verbosity"],
)
)
# Turn into dict
multiplicity_fraction_dict = {
el + 1: multiplicity_fraction_list[el]
for el in range(len(multiplicity_fraction_list))
}
# Set the prev value back
self.population_options["Moe2017_options"]["M_1"] = prev_M1_value_ms
# we don't know what to do next
else:
msg = "Chosen value for the multiplicity fraction function is not known."
raise ValueError(msg)
# To make sure we normalize the dictionary
multiplicity_fraction_dict = normalize_dict(multiplicity_fraction_dict)
self.vb_debug(
"Multiplicity: {} multiplicity_fraction: {}".format(
system_dict["multiplicity"],
multiplicity_fraction_dict[system_dict["multiplicity"]],
),
)
return multiplicity_fraction_dict[system_dict["multiplicity"]]
[docs] def get_moe_di_stefano_dataset(self, options, verbosity=0):
"""
Function to get the default Moe and di Stefano dataset or accept a user input.
Returns a dict containing the (JSON) data.
"""
json_data = None
if "JSON" in options:
# use the JSON data passed in
json_data = options["JSON"]
elif "file" in options:
# use the file passed in, if provided
if not os.path.isfile(options["file"]):
self.vb_error(
"The provided 'file' Moe and de Stefano JSON file does not seem to exist at {}".format(
options["file"]
),
)
raise ValueError
if not options["file"].endswith(".json"):
self.vb_error(
"Provided filename is not a json file",
)
raise ValueError
else:
# Read input data and Clean up the data if there are white spaces around the keys
with open(options["file"], "r", encoding="utf-8") as data_filehandle:
datafile_data = data_filehandle.read()
datafile_data = datafile_data.replace('" ', '"')
datafile_data = datafile_data.replace(' "', '"')
datafile_data = datafile_data.replace(' "', '"')
json_data = json.loads(datafile_data)
if not json_data:
# no JSON data or filename given, use the default 2017 dataset
self.vb_info(
"Using the default Moe and de Stefano 2017 datafile",
)
json_data = copy.deepcopy(moe_di_stefano_2017_data.moe_di_stefano_2017_data)
return json_data